# Wektor wartości zmiennej badanej w populacji
N=25 # liczebność populacji
set.seed(123)
y=round(rnorm(N, mean=4000, sd=800),2)

# Wartość parametru (średnia) w populacji
teta_sr<-mean(y)
teta_sr

# losowanie próby
n=floor(0.2*N)
s <- sample(1:N, size=n, replace=T)
s
proba <- y[s]
proba

# wartość estymatora
est_sr<-mean(proba) # median(proba) # (min(proba)+max(proba))/2 # quantile(proba, probs=0.75)[[1]]
est_sr

# błąd szacunku
u<-est_sr-teta_sr
u

# wszystkie możliwe próby
S=combn(N,n)

# wektory pomocnicze do zapisu wartości błędów etymacji i wartości estymatorów
u_ciag<-rep(NA,dim(S)[2])
est_sr_ciag<-rep(NA,dim(S)[2])

# W pętli for wyznaczamy wartości estymatora dla każdej możliwe próby
for (i in 1:dim(S)[2]) {
	proba<-y[S[,i]]
	est_sr_ciag[i]<- mean(proba) # median(proba) # (min(proba)+max(proba))/2 # quantile(proba, probs=0.75)[[1]]
	u_ciag[i]<-est_sr_ciag[i]-teta_sr
}

# symulacyjna wartość obciązenia estymatora średnej
B_est=mean(u_ciag)
round(B_est,3)

# symulacyjna wartość wariancji estymatora średnej
D2_est <-mean((est_sr_ciag-mean(est_sr_ciag))^2)
D2_est
# średni błąd szacunku
D_est <- sqrt(D2_est)
D_est

# symulacyjna wartość błędu średniokwadratwego estymatora średniej
MSE_est <-mean(((est_sr_ciag-teta_sr))^2)
MSE_est
# pierwiastek z błędu średniokadratowego
RMSE_est <- sqrt(MSE_est)
RMSE_est